
\chapter{蒙特卡洛方法}
我个人相信，人类的博彩业可以追溯到史前文明，每当色(shai三声)子被掷下，金钱就从一个人流向另一个人。
在欧洲的摩纳哥公国，有一个著名的城市叫做蒙特卡洛，每天都在进行大量的赌博。
无数个色子在光洁的桌面上滚动，伴随着男男女女的欢呼与悲叹。
光洁的色子上，镌刻着人类命运的点数，随机地倒向阴暗与光明。
随机性，是那么让人痴迷！
于是聪明绝顶的冯洛伊曼，决定用蒙特卡洛这个城市的名字，命名一类用到随机数的方法。

\section{随机数的产生，乘加同余法}
\subsection{伪随机数}
传统的计算机进行的操作，都是决定性的，如果程序没有写错，给定输入条件，就会得到确定的输出结果。
所以传统的计算机不产生随机性。
于是我们用计算机产生一个无限长的数列，让数列上的数彼此间没有关联，整体均匀分布在某个区间，需要随机数的时候，从中一个一个取用。
这些数叫做“伪随机数”，比较常用的公式是乘加同余法，
\begin{eqnarray}
I_{j+1} = a I_j + c (\rm mod~~ m),
\end{eqnarray}
其中 $a,c,m$ 均为正整数，产生均匀分布的随机整数，当$m$为很大的整数时具有较大的周期。

为了叙述方便，后文中我们不再强调“伪随机数”，而只说“随机数”。
在c++中，cstdlib库中的rand()函数产生[0,RAND\_MAX]区间内的随机数，RAND\_MAX的值是一个很大的整数，与所用的机器有关。
下面的代码输出RAND\_MAX的值，并产生10个随机数。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cstdlib>

int main(){

        int i;
        cout<<"RAND_MAX="<<RAND_MAX<<endl;
        for(i=0;i<10;i++)
                cout<<rand()<<endl;
        return 0;
}
\end{lstlisting}
如果多运行几次上面的程序，你会发现每次输出结果都是一样的，显得不是很“随机”。
这是因为每次运行伪随机数序列的初始值（即“种子”）是一样的，如果想增加一点随机性，可以增加一行，将随机数种子设定为1970年1月1日0点到运行程序的那一刻的总秒数，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cstdlib>

int main(){

        srand(time(NULL)); // 随机数种子：1970年1月1日0点到此刻的总秒数
        int i;
        cout<<"RAND_MAX="<<RAND_MAX<<endl;
        for(i=0;i<10;i++)
                cout<<rand()<<endl;
        return 0;
}
\end{lstlisting}
我们可以用rand()构造离散的随机分布，也可以构造近似的连续随机分布。
在下面的随机点名程序中，我们用rand()对15求余数，得到在$0,1,\cdots,14$中的离散随机数，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<string>
#include<time.h>
#include<vector>

const string name[15]={"邓琦","李雪","刘浩","吕彦萌","王庆轩","王群","王小鹏","王妍","杨浩","杨红艳","张晓璐","赵翔宇","甄洪悦","周琪琪","张莹"};

int main(){

        int i,j,n=-1;
        srand(time(NULL));

        cout<<"你要点多少次名？ n="; cin>>n;
        for(j=0;j<n;j++){
                i = rand() % 15;
                cout<<j+1<<"\t"<<name[i]<<endl;
        }
        return 0;
}
\end{lstlisting}
我们也可以构造$[0,1]$区间内随机分布的“连续”随机数，
\begin{lstlisting}
(double)rand()/RAND_MAX
\end{lstlisting}

\subsection{产生任意概率分布的随机数}
如果你需要一些随机数，分布在$x\in[a,b]$区间内，其概率密度函数为$g(x)$，且有归一化公式
\begin{equation}
\int^b_a g(x) dx = 1,
\end{equation}
那么你可以如下构造
\begin{shaded*}
\begin{itemize}
\item [1] 找到一个值$Max$，使得$\forall x \in [a,b]$，$g(x) \leq Max$。
\item [2] 构造随机数$\xi$，在$[0, Max]$均匀分布
            \lstset{numbers=none}
            \begin{lstlisting}
            (double)rand() / RAND_MAX * Max
            \end{lstlisting}
\item [3] 构造随机数$x$，在$[a,b]$均匀分布
            \lstset{numbers=none}
            \begin{lstlisting}
            (double)rand() / RAND_MAX * (b-a) + a
            \end{lstlisting}
\item [4] 每取一个$x$，也取一个$\xi$；如果$\xi < g(x)$就使用这个$x$；如果$\xi > g(x)$就不使用这个$x$。
\end{itemize}
\end{shaded*}
我们可以把这个算法写成一个函数，供后面调用
\begin{lstlisting}
/*
 * random(...)  产生在 [a,b] 区间内，以 func(x) 为概率密度的随机数，在使用这个函数之前，可以用srand()
函数对随机数种子进行初始化
 * double a
 * double b                     随机数所在区间 [a,b]
 * double (*func)(double x)     随机数的相对概率密度
 * double fmax                  相对概率密度的上界
 */
double random(double a, double b, double (*func)(double x), double fmax){
        
        double x = (double)rand() / RAND_MAX * (b-a) + a; // 在 [a,b] 之间均匀分布
        double f = (double)rand() / RAND_MAX * fmax; // 在 [0, fmax] 之间均匀分布
        while( f > func(x) ){ // x 有 func(x) 的相对概率被接受
                x = (double)rand() / RAND_MAX * (b-a) + a;
                f = (double)rand() / RAND_MAX * fmax;
        }
        return x;
}
\end{lstlisting}
另外我们还可以写一个函数，对随机数在各个区间出现的频率进行统计
\begin{lstlisting}
/*
 * count_random(...) 统计随机数rand出现在 [a,a+h), [a+h, a+2h), ... 中的频率，记录在count中
 */
void count_random(double a, double h, int * count, double rand){
        
        count[ (int)((rand-a)/h) ] ++;
}  
\end{lstlisting}
我们可以产生$[-3,3]$区间内服从标准高斯分布$N(\mu=0, \sigma^2=1)$的随机数
\begin{equation}
f(x) = \frac{1}{\sqrt{2\pi}} exp(-\frac{x^2}{2}),
\end{equation}
并且统计产生的随机数，与概率密度函数相比照，检验我们的代码
\begin{lstlisting}
头文件。。。

double f(double x){//标准高斯分布，或标准正态分布
        return exp(-x*x/2)/sqrt(2*M_PI);
}

int main(){

        int i,N=1000000;
        double a=-3, b=3, fmax = 1/sqrt(2*M_PI);
        double h=0.1, x;
        int n = (b-a)/h;
        int *count = new int [n];
        for(i=0;i<n;i++) count[i] = 0;
        for(i=0;i<N;i++){
                x = gnr_rand(a, b, f, fmax);
                count_random(a, h, count, x);
        }
        ofstream fp("Gauss_rand.txt");
        for(i=0;i<n;i++){
                fp<< a+ (i+0.5)*h <<"\t"<<(double)count[i]/N/h<<endl;
                //第i个区间内的数的频率除以区间长度，得到“频率密度”
        }
        fp.close();
        delete [] count;
        return 0;
}
\end{lstlisting}
得到“频率密度”数据，以后，可以与概率密度函数相比较，gnuplot代码为
\begin{lstlisting}
plot exp(-x*x/2)/sqrt(2*3.1415926), "Gauss_rand.txt" w p
\end{lstlisting}
得到图片\ref{fig:Gauss_rand}，可以看到，我们产生的随机数很好地符合高斯分布。
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{Gauss_rand}
\caption{标准高斯分布，随机数“频率密度”与概率密度的比照}
\label{fig:Gauss_rand}
\end{figure}
如果需要其他分布的随机数，只需修改程序中的f(x)和fmax即可。

\subsection{大数定律与中心极限定理}
随机数$X$服从任意一个分布，其概率密度函数为$f(x)$，期望值为
\begin{equation}
E(X) = \int^\infty_{-\infty} x f(x) dx,
\end{equation}
方差为
\begin{equation}
D(X) = \int^\infty_{-\infty} (x - E(X))^2 f(x) dx = \sigma^2.
\end{equation}
那么，如果每次抽样取$n$个独立服从$f(x)$的随机数，然后取其平均数，则有
\begin{eqnarray}
\lim_{n \rightarrow \infty} P \left\{ | \frac{1}{n} \sum\limits^n_{i=1} X_i - \mu | \leq \epsilon \right\}
\geq \lim_{n \rightarrow \infty} \left[1- \frac{ \sigma^2 }{ n \epsilon^2} \right]
 =1.
\end{eqnarray}
这就是{\bf 大数定律},它的意思很直白：随着抽样次数的增大，抽样平均一定会收敛于期望值。
方差越小，收敛得越快。

{\bf 中心极限定理}则不那么直白，但是可能更神奇：
\begin{eqnarray}
\lim_{n \rightarrow \infty} \frac{ \frac{1}{n}\sum\limits^n_{i=1} X_i - \mu }{ \sigma / \sqrt{n} } \sim N(0,1),
\end{eqnarray}
它的意思是说：如果抽样次数足够大，抽样平均服从正态分布$N(\mu, \sigma^2/n)$。
这就很神奇了，不论$X$服从什么分布，不论$f(x)$是什么样的概率密度分布函数，$X$的抽样平均都服从正态分布。
所以正态分布在大自然中具有特殊的意义，它是分布中的分布。

这两个定理在后文中我们会用到，现在我们不妨用刚学会的伪随机数，稍微检验这两个定理。

我们可以使用$[-3,3]$之间的均匀分布，做30次抽样，取平均，然后统计它的“频率”密度，再与$N(\mu, \sigma^2/30) = N(0,0.1)$做比较，从比较图上便可以一目了然地看到，符合得好不好。
下面的代码实现这个比较
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cstdlib>
#include<cmath>
#include<fstream>

/*
 * random(...)  产生在 [a,b] 区间内，以 func(x) 为概率密度的随机数，在使用这个函数之前，可以用srand()函数对随机数种子进行初始化
 * double a
 * double b                     随机数所在区间 [a,b]
 * double (*func)(double x)     随机数的相对概率密度
 * double fmax                  相对概率密度的上界
 */
double gnr_rand(double a, double b, double (*func)(double x), double fmax){

        double x = (double)rand() / RAND_MAX * (b-a) + a; // 在 [a,b] 之间均匀分布
        double f = (double)rand() / RAND_MAX * fmax; // 在 [0, fmax] 之间均匀分布
        while( f > func(x) ){ // x 有 func(x) 的相对概率被接受
                x = (double)rand() / RAND_MAX * (b-a) + a;
                f = (double)rand() / RAND_MAX * fmax;
        }
        return x;
}

double sample(int n, double a, double b, double (*func)(double x), double fmax){

        double y=0;
        for(int i=0;i<n;i++){
                y += gnr_rand(a, b, func, fmax);
        }
        return y/n;
}

/*
 * count_random(...) 统计随机数rand出现在 [a,a+h), [a+h, a+2h), ... 中的频率，记录在count中
 */
void count_random(double a, double h, int * count, double rand){

        count[ (int)((rand-a)/h) ] ++;
}

double f(double x){//[-3:3]之间均匀分布
        if(x >= -3 && x <= 3) return 1.0/6;
        else return 0;
}

double NormDistr(double mu, double sigma, double x){
        return 1 / sqrt(2*M_PI) / sigma * exp( - (x-mu)*(x-mu)/2/sigma/sigma );
}

int main(){

        int i,N=1000000;
        double a=-3, b=3, fmax = 1.0/6;
        double h=0.01, x;
        int n = (b-a)/h;
        int *count = new int [n];
        for(i=0;i<n;i++) count[i] = 0;
        for(i=0;i<N;i++){
                x = sample(30, a, b, f, fmax);
                count_random(a, h, count, x);
        }
        ofstream fp("test_LLNandCLT.txt");
        for(i=0;i<n;i++){
                fp<< a+ (i+0.5)*h <<"\t"<<(double)count[i]/N/h<<"\t"<< NormDistr(0, sqrt(3.0/30), a+ (i+0.5)*h )<<endl;
        }
        fp.close();
        delete [] count;

        return 0;
}
\end{lstlisting}
得到数据文件以后，作图
\begin{lstlisting}
plot "test_LLNandCLT.txt" u 1:2, "" u 1:3 w l
\end{lstlisting}
得到图片\ref{fig:LLNandCLT}
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{test_LLNandCLT}
\caption{检验大数定律与中心极限定理。用$[-3,3]$之间均匀分布的随机数，进行30次抽样并做平均，根据中心极限定理，这个抽样平均服从$N(0,0.1)$分布。绿色折线为$N(0,0.1)$的概率密度函数，蓝色点为抽样平均值的“频率”密度。}
\label{fig:LLNandCLT}
\end{figure}

\subsection{作业}
1. 用点名程序点一百万次程序，并统计每个人被点到的次数，做一个散点图，并观察是否每个人被点到的概率差不多。

\noindent 
2. 用rand()产生一百万个$[0,1]$区间内的随机数，并统计它们分布于$[0,0.1],\cdots,[0.9,1]$区间内的概率，观察是否都差不多。

\section{随机模拟}
我们介绍几个简单的随机模拟。
\subsection{气体扩散}
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{mc_gas}
\caption{蒙特卡洛模拟：气体扩散的简化模型：任一时刻任一气体分子都有一半可能穿过小孔到另一侧。}
\label{fig:mc_gas}
\end{figure}
有一个封闭腔，如图\ref{fig:mc_gas}，初始时，有1千个气体分子在隔板左侧，此后，任意时刻随机选取一个气体分子穿过小孔到达另一边，那么一段时间之后，我们可以预期，气体分子在两边的分布会是均匀的。
我们写一个代码做这个简单的模拟，练习使用这种随机性，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>

int main(){

        int i, j, n=1000, Nt = 10000;
        int * walker = new int [n];
        for(i=0;i<n;i++){
                walker[i] = -1;// -1/1 for left/right
        }
        int count[2]={1000,0};
        ofstream fp("mc_gas.txt");
        for(j=1;j<=Nt;j++){
                i = rand() % n ;
                walker[i] *= -1;
                if( walker[i] == 1 ){
                        count[1] ++; count[0] --;
                }
                else{
                        count[0] ++; count[1] --;
                }
                fp<<j<<"\t"<<count[0]<<"\t"<<count[1]<<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
得到数据以后，用gnuplot如下作图，
\begin{lstlisting}
plot "mc_gas.txt" u 1:2 w p t "number of molecules on the left", 500 w l notitle
\end{lstlisting}
得到图\ref{fig:mc_gas2}，过一段时间以后，两边气体分子数平衡，此后虽然有些涨落，但基本保持平衡。
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{mc_gas2}
\caption{气体扩散简单模拟}
\label{fig:mc_gas2}
\end{figure}

\subsection{链式反应：原子弹临界尺寸}
U235的链式反应为
\begin{equation}
n + U235 = n + n + else,
\end{equation}
其中$else$表示其他反应产物。
反应后的2个中子如果再与2个U235反应，就会产生4个中子，如此下去，很快就会触发大量反应。
每次反应都会释放MeV量级的反应（生化反应是eV量级），所以大量这样的核反应会产生巨大的能量，迅速产生很高的温度，将材料汽化电离，释放冲击波和大量辐射，这就是1945年在广岛长崎发生的原子弹爆炸。

但是，如果一次反应释放的2个中子与U235再次反应的概率不大，那么就不会触发链式反应，就不会爆炸。
所以如果U235材料的体积特别小，一个U235衰变以后释放的2个中子都飞出了材料，就不会爆炸。
如果U235材料的体积足够大，一次反应产生的2个中子基本都会遇到另2个U235，那么就会爆炸。

所以，原子弹有一个临界尺度，超过尺度就会爆炸，小于尺度就是安全的。

我们可以用随机数模拟这个过程。
我们将长度单位选为中子在材料中的平均自由程，所以在这个单位下，平均自由程为1。
为了简化处理，我们假设中子在材料中的自由程（即两次碰撞之间飞行的距离）为$[0,2]$内均匀分布的随机数，那么自然平均自由程是1。

假设有球形原子弹，半径为R。
那么，我们在原子弹体积内均匀取1000个随机点，认为这些点发生了自发裂变，出现了1000个自由中子。
自由中子的飞行方向随机，朝各个立体角都是相同的概率。
所以$(\theta, \phi)$中$\phi$为$[0,2\pi]$内均匀分布的随机数，而$\theta$要使得
\begin{equation}
d \Omega = \sin \theta d \theta d\phi = d \cos \theta d \phi,
\end{equation}
均匀分布，所以$\cos \theta$是$[-1,1]$内均匀分布的随机数。
那么，自由飞行之后的坐标为
\begin{eqnarray}
x^\prime &=& x + r \sin \theta \cos \phi, \\
y^\prime &=& x + r \sin \theta \sin \phi, \\
z^\prime &=& x + r \cos \theta, \\
\end{eqnarray}
终点若在弹体内部，则有
\begin{equation}
(x^\prime)^2 + (y^\prime)^2 + (z^\prime)^2 \leq R^2,
\end{equation}

我们可以分别计算这1000个中子的自由程；如果落在原子弹内部，就触发一次核反应，增加2个中子；如果落在原子弹外部，就消失不见了，减少1个中子。
最终我们统计中子总数，除以1000，如果这个比例大于1，就意味着下一轮核反应后中子数会继续增大，直到爆炸。
所以比例为1时，对应的就是球形原子弹的临界尺寸，单位是平均自由程。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

#include"library.h"

const int n = 1000;

double ratio( double R, int n ){
        
        int i, count = n;
        double r, costheta, phi, x, y, z;

        srand(time(NULL));
        for(i=0;i<n;i++){

                r = (double)rand()/RAND_MAX * R;
                costheta = (double)rand()/RAND_MAX * 2 -1;
                phi = (double)rand()/RAND_MAX * 2 * M_PI;
                x = r * sqrt( 1 - costheta * costheta ) * cos(phi);
                y = r * sqrt( 1 - costheta * costheta ) * sin(phi);
                z = r * costheta;
                
                r = (double)rand()/RAND_MAX * 2;
                costheta = (double)rand()/RAND_MAX * 2 -1;
                phi = (double)rand()/RAND_MAX * 2 * M_PI;
                x += r * sqrt( 1 - costheta * costheta ) * cos(phi);
                y += r * sqrt( 1 - costheta * costheta ) * sin(phi);
                z += r * costheta;
                
                if( x*x + y*y + z*z <= R*R ){
                        count ++;
                }       
                else{
                        count --;
                }       
        }       
        return (double)count/n;
}       

double f(double R){
        return ratio(R, n) - 1;
}

int main(){
        
        double Rc = bisection(f, 0, 10, 1E-6);
        cout<<"The critical radius R = "<<Rc<<" (unit is mean free path of neutron)"<<endl;
        return 0;
}     
\end{lstlisting}
可以得到，临界半径大约是1.15个平均自由程。

\section{蒙特卡洛积分}
\subsection{置信度与置信区间}
如果一个随机数服从正态分布$N(\mu, \sigma^2)$，那么它的概率密度函数是
\begin{equation}
f(x) = \frac{1}{\sqrt{2\pi}\sigma} exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right).
\end{equation}
那么，随机抽样落在$[\mu-\sigma, \mu+\sigma]$内的概率为
\begin{eqnarray}
F(\sigma) &=& \int^{\mu+\sigma}_{\mu-\sigma} \frac{1}{\sqrt{2\pi}\sigma} exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right) dx \\
&=& \int^{1}_{-1} \frac{1}{\sqrt{2\pi}} exp\left(-\frac{(x)^2}{2} \right) dx \\
&=& 68.268949 \%
\end{eqnarray}
类似地，随机抽样落在平均数左右2个标准差内的概率为$95.449974\%$，落在左右3个标准差内的概率为$99.730020\%$，落在左右4个标准差内的概率为$99.993666\%$。

在历史上，很多物理量的测量和计算给出的是铁板钉钉的区间$[a-\epsilon, a+\epsilon]$，告诉世人这个物理量一定落在这个区间内。
但是有时候我们做不到这一点，尤其是在蒙特卡洛方法中，我们的计算结果$X$可能是随机的，服从分布$N(\mu, \sigma^2)$，可能理论上我们知道$\mu$就是物理量的理论值，但我们不知道它是多少，我们手上只有$n$次抽样的结果$x_i$，根据中心极限定理，我们知道抽样平均$\bar{x}$服从分布$N(\mu, \sigma^2/N)$，那么我们知道，抽样平均落在$[\mu - \sigma/\sqrt{n}, \mu + \sigma/\sqrt{n}]$的概率是68.268949\%，而
\begin{equation}
\mu - \sigma/\sqrt{n} \leq \bar{x} \leq \mu + \sigma/\sqrt{n}
\end{equation}
等价于
\begin{equation}
\bar{x} - \sigma/\sqrt{n} \leq \mu \leq \bar{x} + \sigma/\sqrt{n},
\end{equation}
所以我们估计的结果是：
\begin{shaded*}
{\color{blue}$\mu \in [\bar{x} - \sigma/\sqrt{n}, \bar{x} + \sigma/\sqrt{n}]$的概率是68.268949\%}。
\end{shaded*}

68.268949\%的可能性，似乎不太靠谱，所以我们也可以改说：
\begin{shaded*}
{\color{blue}$\mu \in [\bar{x} - 4\sigma/\sqrt{n}, \bar{x} + 4\sigma/\sqrt{n}]$的概率是99.993666\%}。
\end{shaded*}
99.993666\%看起来挺可靠了，所以我们希望把区间缩小，这个结论就更强大。
这当然可以通过增大$n$（抽样次数）来实现。

这个区间就叫做{\bf 置信区间}，相应的概率就叫做{\bf 置信度}。

\subsection{简单蒙卡积分}
如果我们要求数值积分
\begin{equation}
I_f = \int^b_a f(x) dx,
\end{equation}
我们可以用前面讲到的辛普生积分、高斯积分等方法。

我们也可以用随机数的办法，来估量这个积分。
我们产生$[a,b]$内均匀分布的随机数，做$N$次抽样得到$x_i$，然后计算$f_i \equiv f(x_i)$，再求平均。
可以想象，如果随机数足够多，那么$[a,b]$内几乎所有小区域都被随机数均匀地覆盖到，那么应有
\begin{equation}
\frac{1}{N} \sum_i f_i \approx \frac{1}{b-a}I_f,
\end{equation}
所以有
\begin{equation}
I_f \approx (b-a) \bar{f},
\end{equation}
其中$\bar{f} = \frac{1}{N} \sum_i f_i$。

你也可以把$f_i$看做是一个随机抽样，那么根据中心极限定理，$N$比较大的时候，$\bar{f}$服从$N(E(f), D(f)/N)$正态分布。

$f_i$的期望值是
\begin{equation}
E(f) = \frac{1}{b-a} \int^b_a f(x) dx,
\end{equation}
那么它的方差是
\begin{equation}
D(f) = \frac{1}{b-a} \int^b_a (f(x)-E(f))^2 dx = \frac{1}{b-a} \int^b_a f(x)^2 dx - E(x)^2.
\end{equation}
我们常常用抽样结果得到$D(f)$的近似估计，即
\begin{equation}
D(f) = \frac{1}{b-a} \int^b_a f(x)^2 dx - E(x)^2 \approx \overline{f^2} - \overline{f}^2.
\end{equation}
这样，我们就可以用前面介绍的置信区间、置信度那一套办法，得到对积分的估计：
\begin{shaded*}
{\color{blue} $I_f \in (b-a)[\bar{f} - 4 \sqrt{ \frac{\overline{f^2} - \overline{f}^2}{n} }, \bar{f} + 4 \sqrt{ \frac{\overline{f^2} - \overline{f}^2}{n} }]$的可能性是99.993666\%。}
\end{shaded*}
这样，我们就得到了$I_f$的估计。
不出太大的意外的话，它都落在上面的区间内，我们可以多进行一些抽样，把区间限定得小一点，这个估计就很有用了。

\subsection{多重蒙特卡洛积分}
如果是多重积分
\begin{equation}
I_f = \int_V f dV,
\end{equation}
我们也可以用蒙特卡洛积分的方法，上文中的公式几乎是不变的，
\begin{shaded*}
{\color{blue} $I_f \in V[\overline{f} - 4 \sqrt{ \frac{\overline{f^2} - \overline{f}^2}{n} }, \bar{f} + 4 \sqrt{ \frac{\overline{f^2} - \overline{f}^2}{n} }]$的可能性是99.993666\%。}
\end{shaded*}
这个区间的长度正比于$\sqrt{1/n}$，与积分重数没有关系，这给蒙特卡洛方法带来优势。
这意味着，只要抽样次数一致，积分重数增长时，积分精度不变。
而在传统的积分方法中，要得到同样的精度，取点数按一维情况下的$m$次幂增长，$m$为积分重数。

\subsection{重要抽样}
对于振荡剧烈的被积函数，如果在积分区域内均匀抽样，$\sqrt{ \overline{f^2} - \overline{f}^2 }$ 的值可能很大，导致置信区间很大，也就是精度不高。
我们可以把原来的积分稍作变形
\begin{equation}
I_f = \int_V f dV = \int_V \frac{f}{g} gdV,
\end{equation}
其中$g$是位置的函数。
这可以这么理解：将每个体积微元扩大为原来的$g$倍，将该处的函数值$f$缩小为$f/g$，那么整体的积分不变。
然后，我们可以进行随机抽样，得到$(f/g)_i$的值。
相应地，我们是在$gdV$这样的体积微元中抽样，相当于在$dV$中以概率密度函数$g$进行抽样，即在$V$中有的地方取点多一点，有的地方取点少一点。

那么，相应地，我们得到积分的数值估计为
\begin{shaded*}
{\color{blue} $I_f \in \int_V g dV[\overline{f/g} - 4 \sqrt{ \frac{\overline{(f/g)^2} - \overline{(f/g)}^2}{n} }, \overline{f/g} + 4 \sqrt{ \frac{\overline{(f/g)^2} - \overline{(f/g)}^2}{n} }]$的可能性是99.993666\%。}
\end{shaded*}
如果我们选取$g$，使得$f/g$的方差比较小（实际上，就是$g$尽量接近$f$的形状)，就会得到精度比较高的结果。

\section{马尔科夫(Markov)链}
根据统计系综理论，正则系综（恒温系统）力学量的宏观测量值，等于不同微观状态做如下统计平均得到的结果，
\begin{eqnarray}
\langle A \rangle = \frac{ \int A e^{-\beta E} d \Omega }{ \int e^{-\beta E} d \Omega},
\end{eqnarray}
其中 $E$ 为微观状态对应的体系总能量。
如果我们构造不同的微观状态，并以$e^{-\beta E}$的相对概率进行抽样，那么如果抽样遍布所有可能的围观状态，则有
\begin{eqnarray}
\langle A \rangle = \frac{1}{n} \sum_i A_i,
\end{eqnarray}
我们当然不可能抽取所有的微观状态。
我们能做的是构造一个子空间，其中有一些微观状态，它们出现的相对概率满足正则分布。
那么则有
\begin{eqnarray}
\langle A \rangle \approx \frac{1}{n} \sum_i A_i.
\end{eqnarray}

如果子空间中的微观状态构成一个序列，这个序列就叫做{\bf 马尔科夫（Markov）链}。
从马尔科夫链上抽样并做平均，计算体系的宏观物理量，这样的模拟就叫做{\bf 麦氏模拟(Metropolis)}。

\subsection{细致平衡}
怎样构造马尔科夫链呢？
我们从一个微观状态$X_1$出发，不断通过微小的改变，得到$X_2, X_3, \cdots$。
假定微观状态$X$在链上第$t$个节点出现的概率为$\rho(X,t)$，那么则有
\begin{equation}
\rho(X,t+1) - \rho(X,t) = \sum_{X'}\rho(X',t)T_{X' \rightarrow X} 
- \sum_{X,t}\rho(X,t)T_{X \rightarrow X'},
\end{equation}
其中$T_{X' \rightarrow X}$是链上$X'$通过微小的改变，变为下一个节点$X$的概率，当然有归一化条件
\begin{equation}
\sum_{X} T_{X' \rightarrow X} = 1.
\end{equation}
如果我们控制了$T_{X' \rightarrow X}$这些“跃迁”概率，那么最终的最终，$\rho(X,t)$一定会达到平衡。
比如，如果我们让所有$T_{X' \rightarrow X}$都相等，那么最终我们会得到各种微观状态的均匀分布。

达到平衡以后，应有$\rho(X,t+1) = \rho(X,t)$，则有
\begin{equation}
\sum_{X'}\rho(X',t)T_{X' \rightarrow X} 
= \sum_{X,t}\rho(X,t)T_{X \rightarrow X'}.
\end{equation}
这个方程的解应该有很多个，但我们可以选用一个简单的：
\begin{equation}
\forall X' \neq X,~~ \rho(X',t)T_{X' \rightarrow X} = \rho(X,t)T_{X \rightarrow X'},
\end{equation}
这个就叫做{\bf 细致平衡条件}。

所以在实践中，如果我们要使得Markov链上微观状态的相对概率为$\rho(X,t)$，那么我们将它带入上式，决定$T_{X \rightarrow X'}$，这样来构Markov链。
经过一段演化之后，Markov链上的分布就会收敛于$\rho(X,t)$。

对于正则系综，由于
\begin{equation}
\frac{\rho(X,t)}{\rho(X',t)} = \frac{e^{-\beta E_X}}{e^{-\beta E_{X'}}}
= e^{-\beta(E_X - E_{X'})},
\end{equation}
所以我们需要构造
\begin{equation}
\frac{ T_{X' \rightarrow X} }{T_{X \rightarrow X‘}} =\frac{\rho(X,t)}{\rho(X',t)} = e^{-\beta(E_X - E_{X'})},
\end{equation}
一种常用的取法是
\begin{eqnarray}
T_{X \rightarrow X'} = \left\{
\begin{aligned}
&e^{ - \beta (E_{X'} - E_{X}) }, &E_{X'} > E_{X};\\
&1, ~~~~~~~~~~~~~~~~~~~~~~~~&E_{X'} < E_{X}.
\end{aligned}
\right.
\end{eqnarray}

\subsection{Ising模型}
一维、二维、三维Ising模型的模拟思路是完全一样的，一维没有相变，三维计算量较大，所以我们展示二维Ising模型的模拟。
二维Ising模型中，原子呈二维晶格，简单起见我们可以设为正方晶格，体系哈密顿量为：
\begin{equation}
H = -J \sum_{\{i,j\}} S_i S_j - \mu B\sum_i S_i,
\end{equation}
其中，$\{i,j\}$表示相邻的两个原子的序号，$S_i$表示第$i$个原子的自旋第3分量，只能取$-1,1$两个值，分别对应自旋向上、自旋向下，$J$表示交换相互作用强度，为正时对应铁磁物质，为负时对应反铁磁物质，$\mu$是原子磁矩，$B$是外加的$z$轴方向匀强磁场。

如果材料属于铁磁物质，则$J$为正，相邻原子取向一致时，体系能量更低，所以整体趋向于所有原子自旋取向一致的情况，这就对应宏观的磁化现象；但是如果温度比较高，非磁化的微观状态也会比较重要，所以随着温度升高，磁性逐渐减弱；在一个特定的温度附近，磁性突然消失，这就是铁磁-非铁磁相变。

如果$J$为负，相邻原子取向相反时，体系能量更低，所以整体趋向于所有院子无规则取向的情况，所以这个材料是反铁磁物质。

如果加外磁场，即$B$不为零时，静磁力使得原子取向趋于某个方向，所以外磁场会影响相变点。

在计算机中，我们用有限大的二维网格模拟Ising模型，Markov链算法描述如下：
\begin{shaded*}
\begin{itemize}
\item[1] 定义$L\times L$的原子点阵，即$n=L^2$个原子呈正方形排布。这些原子的自旋只能朝上或者朝下。
\item[2] 使用周期性边界条件，点阵右边界原子的右邻是左边界原子，上边界原子的上邻是下边界原子。这样，各个原子的地位完全等同。
\item[3] 确定初始状态（最终结果对初始状态不依赖），确定体系温度，相互作用强度。
\item[4] 从初始状态出发，每次随机选一个原子，以如下概率旋转它的自旋
\begin{eqnarray}
T_{X \rightarrow X'} = \left\{
\begin{aligned}
&e^{ - \beta (E_{X'} - E_{X}) }, &E_{X'} > E_{X};\\
&1, ~~~~~~~~~~~~~~~~~~~~~~~~&E_{X'} < E_{X}.
\end{aligned}
\right.
\end{eqnarray}
其中，$X,X'$表示旋转前后整个体系的两个不同微观状态。
\item[5] 经过足够多的旋转以后，比如，Markov链的前1/5过去以后，认为Markov链已经达到平衡，开始取点。
\item[6] 对取出的每个Markov点$i$，即一个微观状态，计算它的物理量$X_i$。
\item[7] 做平均$\frac{1}{N}\sum^N_i X_i$，认为它是系综平均$\langle X \rangle$的近似值，作为结果给出。 
\end{itemize}
\end{shaded*}

我们可以写一个代码，实现上面的算法，同时，将Markov链上选择一些等间距的微观状态，并可视化，制作动图。
下面的代码一边模拟Ising模型，一边选取微观状态，并输出相应的gnuplot指令，生成gnuplot动图脚本
\begin{lstlisting}
/*
 * 2-dimensional Ising model, with Markov chain in Monte Carlo method
 *      H = - J \sum_{ij} s_i s_j - muB\sum_i s_i,
 * where J is the interaction strength between two atoms, negative/positive for ferromagnetic/anti-ferromagnetic material, k is the Boltzman constant, T is the temperature, we set J as unit.
 */

#include<iostream>
#include<stdlib.h>
#include<time.h>
#include<fstream>
#include<cmath>
using namespace std;

int main(){

        long length, i,j, L=4, num=16, S, S_neighbour, count=0;
        double J=1, kT, muB, E, dE, P, S_ave=0, E_ave=0, E2_ave=0;
        int *spin=new int [L*L];//spin of each atom, which can only be 1 or -1
        bool flag_flip;

        for(i=0;i<L*L;i++) //initially the lattice is magnetizied.
                spin[i]=1;

        cout<<"======================================================================"<<endl;
        cout<<"\t\t2D Ising Model for Demonstration\t"<<endl;
        cout<<"----------------------------------------------------------------------"<<endl;

        srand((unsigned)time(0));//initialize rand() with a random seed, according to the real time when this code runs.
        cout<<"\tHamiltonian:\tH= -J \\sum_{ij} S_i S_j - muB\\sum_i S_i"<<endl;
        cout<<"\tParameter:\tlength="; cin>>length;
        cout<<"\t\t\tJ="<<J<<endl;
        cout<<"\t\t\tkT=";  cin>>kT;
        cout<<"\t\t\tmuB=";  cin>>muB;
        cout<<"\t--------------------------------------------------------------"<<endl;

        S=num;//initial value of total spin
        E=-J*2*num;//initial value of total energy of the system

        ofstream fp("IsingModel_test.gnu");
        fp<<"set term gif animate delay 10\n";
        fp<<"set output \"IsingModel.gif\"\n";
        fp<<"unset xtics\n";
        fp<<"unset ytics\n";
        fp<<"unset border\n";
        fp<<"set title \"kT = "<<kT<<"J, n=0, <s>=0, <E>=0\n";
        fp<<"set dummy i,j\n";
        fp<<"do for[i=0:3]{\n";
        fp<<"\t do for[j=0:3]{\n";
        fp<<"        set arrow from j+0.5,i+0.2 to j+0.5,i+0.8 size 3,0.8 filled lw 5 lc rgbcolor \"red\"\n";
        fp<<"\t }\n";
        fp<<"}\n";
        fp<<"plot \"rectangle.txt\" with line lc rgbcolor \"black\" notitle\n";

        for(j=0;j<length;j++){

                if( j% (length/10) ==0 )cout<<"Finished "<<100.0*j/length<<"%.\n";

                i=rand()%num;
                /*
                cout<<"i = "<<i<<endl;
                cout<<"up = "<< (i>11 ? i-12 : i+4) <<endl;
                cout<<"down = "<< (i<4? i+12 : i-4 ) <<endl;
                cout<<"left = "<< (i%4 ==0 ? i+3 : i-1) <<endl;
                cout<<"right = "<< ((i+1)%4==0 ? i-3 : i+1) <<endl;
                */
                S_neighbour = spin[ (i>11 ? i-12 : i+4) ] + spin[ (i<4? i+12 : i-4) ]
                        + spin[ (i%4 ==0 ? i+3 : i-1) ] + spin[ ((i+1)%4==0 ? i-3 : i+1) ];

                dE= 2 * J * S_neighbour * spin[i] + 2 * muB * spin[i];

                flag_flip=false;
                if(dE<0){
                       flag_flip=true;
                }
                else{
                        if( (double)rand()/RAND_MAX < exp(-dE/kT) ){// rand()/RAND_MAX uniformly distributes in [0,1], so it has possibility P to be less than P
                                flag_flip=true;
                        }
                }

                if(flag_flip){
                        E += dE;
                        spin[i] *= -1;//flip the spin
                        S += 2*spin[i];// total spin and total energy changes
                }

                if( fabs((double)S/16) > 1 )cout<<"s="<<(double)S/16<<endl;

                if(j>0.1*length){
                        E_ave = (double)(count)/(count+1) * E_ave + 1.0/(count+1) * E;
                        E2_ave = (double)(count)/(count+1) * E2_ave + 1.0/(count+1) * E * E;
                        S_ave = ((double)count * S_ave + (double)S )/(count+1);
                        count ++;
                        if(j% (length/100)==0){
                                fp<<"set title sprintf(\"kT="<<kT<<"J, n=%d, <s>=%f, <E> = %f\","<<count<<","<<(double)S_ave/16<<","<<E_ave<<")\n";
                                fp<<"unset arrow "<<(i+1)<<"\n";
                                if(spin[i]>0)
                                        fp<<"set arrow from "<<(i%4)+0.5<<","<<i/4+0.2<<" to "<<(i%4)+0.5<<","<<i/4+0.8<<" size 3,0.8 filled lw 5 lc rgbcolor \"red\"\n";
                                else
                                        fp<<"set arrow from "<<(i%4)+0.5<<","<<i/4+0.8<<" to "<<(i%4)+0.5<<","<<i/4+0.2<<" size 3,0.8 filled lw 5 lc rgbcolor \"blue\"\n";
                                fp<<"replot\n";
                        }
                }
                        /* for test only
                        for(int p=0;p<L;p++){
                                for(int q=0;q<L;q++){
                                        if(spin[ p*L + q]==1)cout<<"+\t";
                                        else cout<<"-\t";
                                }
                                cout<<endl;
                        }
                        */
                //fp<<j<<"\t"<<E_ave<<"\t"<<E2_ave<<"\t"<<S_ave<<endl;
        }
        fp<<"set output\n";
        fp.close();
        cout<<"Results:"<<endl;
        cout<<"\t <E> = "<<E_ave<<endl;
        cout<<"\t <E2> = "<<E2_ave<<endl;
        cout<<"\t <s> = "<<(double)S_ave/16<<endl;

        //double C=1/kT/kT*( E2_sample/num_sample - (E_sample*E_sample/num_sample/num_sample) );
        return 0;
}             
\end{lstlisting}
得到gnuplot脚本以后，可以制作动图，在课堂上演示。

\section{变分蒙特卡洛}
\subsection{变分法}
在初等量子力学中，我们学过变分法，用来优化体系的一个近似本征态，下面简单回顾这一理论。
若体系的哈密顿量 $\hat{H}$ 已知，存在本征值 $\epsilon_i$ 与 $\phi_i$，使得
\begin{eqnarray}
\hat{H} \phi_i = \epsilon_i \phi_i,~~ i = 1, 2, 3, \cdots.
\end{eqnarray}
如果问题复杂，无法严格解出 $\epsilon_i, \phi_i$，那么能得到基态能量和基态波函数的近似解，就已经很满足了。

于是我们构造一个波函数 $\Psi$。无论怎样构造，它一定可以表示成 $\phi_i$ 的线性展开（因为 $\hat{H}$ 的所有本征态构成完备基），
\begin{eqnarray}
\Psi = a_1 \phi_1 + a_2 \phi_2 + \cdots,
\end{eqnarray}
虽然我们不知道 $\phi_1, \phi_2, \cdots $ 长什么样子，$\Psi$ 是它们的线性展开，这一点总是没错的。

另外，
\begin{eqnarray}
\langle E \rangle &=& \frac{ \langle \Psi | \hat{H} | \Psi \rangle }{ \langle \Psi | \Psi \rangle }
= \frac{ \sum\limits^n_{i,j=1} a^*_i a_j \langle \phi_i | \hat{H} | \phi_j \rangle }{ \sum_i |a_i|^2 }
= \frac{ \sum\limits^n_{i,j=1} a^*_i a_j \epsilon_j \delta_{ij} }{ \sum_i |a_i|^2 }
\nonumber \\
&=& \frac{ \sum\limits^n_{i=1} |a_i|^2 \epsilon_i }{ \sum_i |a_i|^2 }
\geq \frac{ \sum\limits^n_{i=1} |a_i|^2 \epsilon_1 }{ \sum_i |a_i|^2 }
= \epsilon_1
\end{eqnarray}
当且仅当 $a_2 = a_3 = \cdots =0$ 时，上式中的 $\geq$ 符号取等号。
也就是说，对于任意不为0的试探波函数，$ \langle E \rangle$ 一定不小于体系基态能量。
$| \Psi \rangle$ 越接近体系基态波函数，$ \langle E \rangle$ 就会越接近体系基态能量。
$| \Psi \rangle$ 往往有一个带参数的拟设形式，在此形式下，$ \langle E \rangle$ 往往是可以计算的，所以我们要做的，就是调节 $\Psi$（即调节其中参数），使得 $ \langle E \rangle$ 尽量小，取极小值的时候，就把 $\Psi$ 看作 $\phi_1$ 的近似，把 $ \langle E \rangle$ 看作是 $\epsilon_1$ 的近似。
这就是变分法。

\subsection{变分蒙特卡洛(VMC)}
在给定参数$\alpha$下，我们可以直接做积分，然后计算
\begin{equation}
\langle E \rangle_\alpha = \frac{ \langle \Psi_\alpha | \hat{H} | \Psi_\alpha \rangle }{ \langle \Psi_\alpha | \Psi_\alpha \rangle },
\end{equation}
然后改变$\alpha$，使得上式达到极小值。
$\Psi_\alpha(\vec{r})$在$\vec{r}$空间上往往是变化很大的，在大部分区域接近于0，可以使用重要抽样。

我们也可以用Markov链的思路，换一种方法求取$\langle E \rangle_\alpha$。

首先，在$\vec{r}$区间内随机选择多个随机行走者（random walker）。
对于每一个 random walker，计算
\begin{equation}
E_L (\vec{r})= H \Psi (\vec{r}) / \Psi(\vec{r}),
\end{equation}
这个量叫做局域能量，我们可以写出
\begin{equation}
\langle E \rangle = \frac{ \int \Psi(\vec{r})^2 E_L(\vec{r}) d\vec{r} }{ \int \Psi(\vec{r})^2  d\vec{r} },
\end{equation}
所以，相当于以$\Psi(\vec{r})^2$为相对概率，对$E_L(\vec{r})$进行积分。
所以我们可以构造相对概率为$\Psi(\vec{r})^2$的马尔科夫链，对$E_L(\vec{r})$进行抽样，得到$\langle E \rangle$。

\begin{shaded*}
\begin{itemize}
\item[1] 在$\vec{r}$空间均匀选取大量random walker。
\item[2] 每个 random walker 计算
\begin{equation}
\langle E \rangle_\alpha = \frac{ \langle \Psi_\alpha | \hat{H} | \Psi_\alpha \rangle }{ \langle \Psi_\alpha | \Psi_\alpha \rangle },
\end{equation}
\item[3] 每个 random walker 产生一个随机行走步长，比如$x,y,z$都有一般概率加上$\delta$，一般概率减去$\delta$；也可以做得更均匀一点，让步长为朝随机立体角行走$\delta$。
然后计算终点的波函数值　$\Psi(\vec{r}^\prime)$，并计算比例
\begin{equation}
p = \Psi(\vec{r}^\prime)^2 / \Psi(\vec{r})^2,
\end{equation}
\item[4] 根据如下概率，决定是否接受这次随机行走
\begin{eqnarray}
T_{\vec{r} \rightarrow \vec{r}^\prime} = \left\{
\begin{aligned}
&p, &p<1;\\
&1, &p>1.
\end{aligned}
\right.
\end{eqnarray}
\item[5]根据细致平衡条件，最终random walker出现在各处的概率与该处波函数平方成正比。
所以我们做$E_L$的抽样即可，
\begin{equation}
\langle E \rangle_\alpha \approx \overline{E_L}.
\end{equation}
\item[6]变化$\alpha$的值，重复上述过程，做出$\overline{E_L} \sim \alpha$的曲线，找到极小值点，即完成近似求解。
\end{itemize}
\end{shaded*}

\paragraph{一维谐振子}
没有量纲的一维谐振子薛定谔方程为
\begin{equation}
H \psi(x) = (-\frac{1}{2} \frac{d^2}{dx^2} + \frac{1}{2} x^2) \psi(x),
\end{equation}
基态的严格解为
\begin{equation}
\psi_{g.s.} (x) = A e^{-x^2/2},
\end{equation}
其中$A$为归一化因子。
基态能量为$e_{g.s.} = 1/2$。

为了演示VMC的思路，我们选择拟设形式
\begin{equation}
\psi_\alpha (x) = e^{-\alpha x^2},
\end{equation}
则局域能量为
\begin{equation}
E_L(x) = H \psi(x) / \psi(x) = \alpha + (\frac{1}{2} - 2\alpha^2) x^2,
\end{equation}
利用这个公式，如下代码实现VMC，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>

double psi(double alpha, double x){
        return exp(-alpha*x*x);
}
double EL(double alpha, double x){
        return alpha + (0.5-2*alpha*alpha)*x*x;
}

double vmc1dho_aveEL(double alpha, int n, double h, int N){

        srand(time(NULL));
        int i,j,count;
        double a=-5, b=5, x, d, y, z, p;
        y=0;
        for(i=0;i<n;i++){
                x = (double)rand()/RAND_MAX * (b-a) + a;//初始位置
                z=0;
                count=0;
                for(j=0;j<N;j++){
                        if( j > 0.1 *N ){
                                z += EL(alpha, x);
                                count ++;
                        }
                        d = ((rand() % 2)*2-1)*h;
                        p = pow( psi(alpha, x+d), 2) / pow( psi(alpha,x), 2);
                        if(p>1) x += d;
                        else{
                                if( (double)rand()/RAND_MAX < p ) x+=d;
                        }
                }
                z /= count;
                y += z;
        }
        return y/n;
}

int main(){

        ofstream fp("vmc1dho_aveEL.txt");
        for(double alpha=0.4; alpha <=0.6; alpha += 0.01){
                fp << alpha << "\t" << vmc1dho_aveEL(alpha, 100, 0.2, 10000) <<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
代码运行以后，通过以下gnuplot代码可以生成图片
\begin{lstlisting}
set term png
set output "vmc1dho.png"
set xlabel "{/Symbol a}"
set ylabel "<E>"
plot "vmc1dho_aveEL.txt" u 1:2 w l lw 2 t "VMC: 1d Harmonic Oscillator"
\end{lstlisting}
得到图片\ref{fig:vmc1dho}。
可以看到，最低点对应着$\alpha=1/2$，即体系基态波函数，相应的$\langle E \rangle = 1/2$对应着基态能量$1/2$。
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{vmc1dho}
\label{fig:vmc1dho}
\end{figure}

\paragraph{氢原子}
没有量纲的氢原子薛定谔方程为
\begin{equation}
(-\frac{1}{2}\nabla^2 - \frac{1}{r} ) \psi(\vec{r}) = E \psi(\vec{r}),
\end{equation}
分离变量以后，设
\begin{equation}
\psi(\vec{r}) = ( u(r)/r ) Y_{lm}(\theta, \phi),
\end{equation}
$s$波方程转化为一维径向方程
\begin{equation}
(-\frac{1}{2} \frac{d^2}{dr^2} - \frac{1}{r} ) u(r) = E u(r),
\end{equation}
基态的严格解为
\begin{equation}
u(r) = A r e^{-r},
\end{equation}
对应着
\begin{equation}
\psi_{g.s.} (\vec{r}) = A e^{-r}.
\end{equation}
其中$A$为归一化因子。
基态能量为$e_{g.s.} = -1/2$。

为了演示VMC的思路，我们选择拟设形式
\begin{equation}
u(r) = r e^{-\alpha r},
\end{equation}
没有归一化，但是不影响我们后面的计算，因为我们真正使用的是波函数在不同地方的相对大小。

局域能量为
\begin{equation}
E_L(r) = (-\frac{1}{2} \frac{d^2}{dr^2} - \frac{1}{r} )　u(r) / u(r) = -\frac{1}{2}\alpha^2 + (\alpha-1)/r,
\end{equation}
利用这个公式，如下代码实现VMC，

\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>
#include"library.h"
double u(double alpha, double r){
        return r*exp(-alpha*r);
}
double EL(double alpha, double r){
        return -0.5*alpha*alpha + (alpha-1)/r;
}
int main(){

        ofstream fp("vmcHydrogen_aveEL.txt");
        for(double alpha=0.8; alpha <=1.2; alpha += 0.01){
                fp << alpha << "\t" << vmc1d_aveEL(alpha, 100, u, EL, 0, 5, 0.2, 10000) <<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}

得到下图。
可以看到，最低点对应着$\alpha=1$，即体系基态波函数为$e^{-r}$，相应的$\langle E \rangle = -1/2$对应着基态能量$-1/2$。
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{vmcHydrogen}
\caption{变分蒙特卡洛：氢原子拟设波函数$\psi_\alpha(\vec{r}) = e^{-\alpha r}$，在$\alpha$取$0.8\sim 1.2$时，得到的$<E> = <\psi_\alpha|H|\psi_\alpha>/<\psi_\alpha | \psi_\alpha>$值。}
\label{fig:vmcHydrogen}
\end{figure}

